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ABSTRACT 



We present a new numerical code, X-ECHO, for general relativistic magnetohydrodynamics (GRMHD) in dynamical spacetimes. 
This is aimed at studying astrophysical situations where strong gravity and magnetic fields are both supposed to play an important 
role, such as for the evolution of magnetized neutron stars or for the gravitational collapse of the magnetized rotating cores of massive 
stars, which is the astrophysical scenario believed to eventually lead to (long) GRB events. The code is based on the extension of 
the Eulerian conservative high-order (ECHO) scheme [Del Zanna et ah, A&A 473, 11 (2007)] for GRMHD, here coupled to a novel 
solver for the Einstein equations in the extended conformally flat condition (XCFC). We solve the equations in the 3+1 formalism, 
assuming axisymmetry and adopting spherical coordinates for the conformal background metric. The GRMHD conservation laws 
are solved by means of shock-capturing methods within a finite-difference discretization, whereas, on the same numerical grid, the 
Einstein elliptic equations are treated by resorting to spherical harmonics decomposition and solved, for each harmonic, by inverting 
band diagonal matrices. As a side product, we build and make available to the community a code to produce GRMHD axisymmetric 
equilibria for polytropic relativistic stars in the presence of differential rotation and a purely toroidal magnetic field. This uses the 
same XCFC metric solver of the main code and has been named XNS. Both XNS and the full X-ECHO codes are validated through 
several tests of astrophysical interest. 
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1. Introduction 

The most spectacular phenomena in high-energy Astrophysics, 
like those associated to active galactic nuclei (AGNs), galactic 
X-ray binary systems, or gamma-ray bursts (GRBs), typically 
involve the presence of rotating compact objects and magnetic 
fields. In some cases, such as the merging of binary systems 
(formed by either neutron stars, NSs, or black holes, BHs) or the 
collapse of rotating cores of massive stars towards a NS or BH, 
the interplay between matter, electromagnetic fields, and gravity 
is so strong that the MHD equations governing the fluid motions 
must be solved self-consistently with the Einstein equations for 
the spacetime metric. Even for less violent phenomena, as the 
oscillations of neutron stars (Font et al. 2002), a self consis- 
tent solution of the fluid equations together with the spacetime 
evolution is essential to properly estimates the frequencies of 
the eigenmodes. In the case of a binary NS merger, which is a 
possible mechanism to account for short GRBs, a strong mag- 
netic field could be produced due to the induced shear (Price 
& Rosswog 2006). Long GRBs are instead associated to super- 
nova events and to the core collapse of massive stars (Woosley 
& Bloom 2006), leading to the subsequent formation of a ro- 
tating and strongly magnetized compact object. The mainstream 
collapsar model (Woosley 1993) implies the rapid formation of 
a maximally rotating Kerr BH at the center, accreting material 
from a torus and likely to loose energy through the Blandford- 
Znajeck mechanism (Barkov & Komissarov 2008). However, a 
promising alternative for the GRB central engine involves the 
presence of a millisecond magnetar with B > 10 15 G (Usov 
1992). The origin of such enormous fields for magnetars is prob- 



ably due to efficient dynamo action during the neutrino cooling 
phase in the hot, deleptonizing proto-NS (Duncan & Thompson 
1992). On the observational side, magnetars are the accepted ex- 
planation for anomalous X-ray pulsars and soft gamma-ray re- 
peaters (Kouveliotou et al. 1998). 

On the computational side, the last decade has witnessed 
a very rapid evolution in the construction of shock-capturing 
codes for general relativistic MHD (GRMHD) in both static 
and dynamical spacetimes, with a wealth of astrophysical ap- 
plications to the situations outlined above (Font 2008). In the 
present paper we describe a novel code for GRMHD in dynam- 
ical spacetimes, named X-ECHO, aimed at studying the evolu- 
tion of magnetized relativistic stars and the gravitational collapse 
of the magnetized rotating cores of massive stars. X-ECHO is 
built on top of the Eulerian conservative high-order code (Del 
Zanna et al. 2007) for GRMHD in a given and stationary back- 
ground metric (Cowling approximation), which in turn upgraded 
the previous version for a Minkowskian spacetime (Del Zanna & 
Bucciantini 2002; Del Zanna et al. 2003). ECHO relies on robust 
shock-capturing methods within a finite-difference discretization 
scheme (two-wave Riemann solvers and limited high order re- 
construction routines), with a staggered constrained transport 
method to preserve the divergence-free condition for the mag- 
netic field (to machine accuracy for second order of spatial accu- 
racy), as proposed by Londrillo & Del Zanna (2000, 2004). The 
ECHO code has been already successfully applied to a variety 
of astrophysical situations involving magnetized plasmas around 
compact objects, like the dynamics and non-thermal emission of 
pulsar wind nebulae (Bucciantini et al. 2003; Del Zanna et al. 
2004; Bucciantini et al. 2004, 2005a,b; Del Zanna et al. 2006; 
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Volpi et al. 2008), emission of relativistic MHD winds from ro- 
tating NSs (Bucciantini et al. 2006), magnetar winds producing 
long GRB jets escaping the stellar progenitor (Bucciantini et al. 
2008, 2009), and post-merger accreting disks around Kerr BHs 
(Zanotti et al. 2010). In spite the code being fully 3D, due to the 
nature of the sources, invariably a plasma surrounding a central 
compact object, all the above applications have been performed 
in 2D axisymmetric spacetimes using spherical-type coordinates 
(either in Minkowski, Schwarzschild, or Kerr metric). The X- 
ECHO version presented and tested here shares the same philos- 
ophy, and, in view of the future applications mentioned above, 
only the axisymmetric case will be considered. 

The Einstein and GRMHD equations in X-ECHO are written 
by fully exploiting the so-called 3 + 1 formalism (like in ECHO), 
in which the original equations are split in their temporal and 
spatial components. The 3 + 1 formalism is nowadays adopted in 
basically all numerical schemes for general relativity (Alcubierre 
2008; Baumgarte & Shapiro 2010), where the system of Einstein 
equations is treated like a Cauchy problem with some initial data 
to be evolved in time through hyperbolic equations. However, 
like for the solenoidal condition for the magnetic field, non- 
evolutionary constraints must be preserved in the numerical evo- 
lution and computational methods for modern codes are divided 
into two main classes: 1) free-evolution schemes, mainly based 
on hyperbolic equations alone, where this problem is alleviated 
by appropriate reformulations of the equations (BSSN: Shibata 
& Nakamura 1995; Baumgarte & Shapiro 1999), eventually with 
the addition of propagating modes and damping terms (Z4: Bona 
et al. 2003; Bernuzzi & Hilditch 2010); 2) fully constrained 
schemes, where the constraints are enforced at each timestep 
through the solution of elliptic equations (Bonazzola et al. 2004), 
a more robust but computationally demanding option, since el- 
liptic solvers are notoriously difficult to parallelize. Most of the 
state-of-the-art 3D codes for GRMHD in dynamical spacetimes 
are based on free-evolution schemes in Cartesian coordinates 
(Duez et al. 2005; Shibata & Sekiguchi 2005; Anderson et al. 
2006; Giacomazzo & Rezzolla 2007; Montero et al. 2008; Farris 
et al. 2008), and have been used, in the case of magnetized plas- 
mas, for gravitational collapse (Duez et al. 2006a; Shibata et al. 
2006a,b; Stephens et al. 2007, 2008), evolution of NSs (Duez 
et al. 2006b; Kiuchi et al. 2008; Liebling et al. 2010), binary 
NS mergers (Anderson et al. 2008; Liu et al. 2008; Giacomazzo 
et al. 2009, 2010), and accreting tori around Kerr BHs (Montero 
et al. 2010). 

Provided the emission of gravitational waves is not of pri- 
mary interest, a good option in the class of fully constrained 
scheme is represented by the conformally flat condition (CFC) 
schemes (e.g. Wilson & Mathews 2003; Isenberg 2008), an ap- 
proximation often employed for the study of gravitational col- 
lapse or NS stability and evolution (Dimmelmeier et al. 2002; 
Saijo 2004; Dimmelmeier et al. 2006; Cerda-Duran et al. 2008; 
Abdikamalov et al. 2009). CFC is typically associated to ax- 
isymmetric configurations and spherical coordinates, and it is 
exact in the spherically symmetric case. Deviations from full 
GR solutions in the axisymmetric case for this kind of appli- 
cations have already been shown to be negligible (Shibata & 
Sekiguchi 2004; Ott et al. 2007), though the nonlinear equations 
may show serious uniqueness problems for highly compact NSs 
or nascent BHs (see Cordero-Carrion et al. 2009, and references 
therein). Moreover, CFC requires the solution of all the ellip- 
tic equations for the metric terms at the same time, usually by 
means of iteration of Poisson solvers, together with the inversion 
of conservative to primitive fluid/MHD variables, another itera- 
tive numerical process. All these difficulties have been resolved 



recently by the extended conformally flat condition (XCFC) for- 
mulation (Cordero-Carrion et al. 2009): all the elliptic equations 
(now eight rather than five) are hierarchically decoupled and lo- 
cal uniqueness is ensured (see also Saijo 2004), thus this will be 
our choice for X-ECHO. 

In the present work we propose and test a new numerical 
solver based on XCFC for an axisymmetric spacetime in confor- 
mally flat spherical-like coordinates. We employ the same nu- 
merical grid used for the evolution of the fluid and magnetic 
quantities through the ECHO scheme, and the Poisson-like equa- 
tions are solved through a hybrid method based on spherical har- 
monics decomposition and direct inversion of band diagonal ma- 
trices, resulting from a second order finite difference discretiza- 
tion of the radial equations, for each harmonic. As a side prod- 
uct, we build and make available to the community a numerical 
code based on XCFC to produce self-consistent GRMHD ax- 
isymmetric equilibria for polytropic relativistic stars in the pres- 
ence of differential rotation and toroidal magnetic fields, here 
named XNS. Both XNS and the full X-ECHO codes are vali- 
dated through several tests of astrophysical interest, including 
accuracy checks in the initial data for various NS equilibrium 
configurations, accuracy in finding the frequencies of their nor- 
mal modes of oscillations, an evolutionary test of migration of 
NS unstable equilibria to stable branches, a test on the stability 
of a differentially rotating, magnetized NS with a toroidal field, 
ID and 2D collapse of an unstable NS toward a BH, and a toy 
collapse of a differentially rotating NS with poloidal fields, as a 
first step towards more realistic magneto-rotational core collapse 
simulations. 

The paper is structured as follows. In Sect. 2 we introduce 
and review the Einstein equations in the 3 + 1 formalism, first 
in their general form and then specialized in the CFC approx- 
imation, and we review the GRMHD equations. In Sect. 3 we 
discuss the new numerical XCFC solver assuming axisymme- 
try and spherical coordinates for the conformal flat 3-metric, 
whereas the description of our novel XNS code for NS initial 
data can be found in Sect. 4. Numerical validation and testing of 
various cases of NS equilibria, oscillations and collapse are pre- 
sented in Sect. 5, while Sect. 6 is devoted to the conclusions. In 
the following we assume a signature {-, +, +, +} for the space- 
time metric and we use Greek letters [i,v,A,... (running from 
to 3) for 4D spacetime tensor components, while Latin letters 

1, j, k,... (running from 1 to 3) will be employed for 3D spatial 
tensor components. Moreover, we set c = G = Mq = 1 and we 
absorb the V47T factors in the definition of the electromagnetic 
quantities. 

2. Basic equations in the 3 + 1 formalism 

In the present section we present the 3 + 1 formalism for 
Einstein equations. Further details can be found in recent books 
and reviews on 3 + 1 numerical relativity (Gourgoulhon 2007; 
Alcubierre 2008; Baumgarte & Shapiro 2010). We briefly dis- 
cuss the constrained evolution schemes, focussing on the elliptic 
CFC and XCFC solvers, and finally, in Sect. 2.3 we review the 
GRMHD equations in 3 + 1 conservative form, as implemented 
in the original ECHO scheme (Del Zanna et al. 2007). 

2.1. The Einstein equations for conformal flatness 

The field equations of general relativity expressed in the 4D fully 
covariant form 

Guv = SnT^y, (1) 
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where G^ v is the Einstein tensor containing the derivatives of 
the metric tensor g^ and is the matter and/or electromag- 
netic energy-momentum tensor, are not appropriate for numer- 
ical computations since time and space are treated on an equal 
footing, whereas one would like to cast them in the form of an 
initial value (or Cauchy) problem and evolve them in time. The 
most widely used approach to reach this goal is based on the so- 
called 3 + 1 formalism, in which the generic spacetime (At, g^y) 
is split into space-like hyper-surfaces Z t . If n 1 * indicates the time- 
like unit normal to S, (also known as the velocity of the Eulerian 
observer, n^rf = -1), the induced three-metric on each hyper- 
surface and the related extrinsic curvature can be defined respec- 
tively as 

7„v -=gny + ri/jlly, (2) 

V := -y^nv, (3) 

where is the covariant derivative with respect to g^y (so that 
^Agfiv = 0). In general, any 4-vector or tensor can be decom- 
posed into normal and spatial components by contracting with 
-n M or with the projector jy, respectively. In particular, both y MV 
and K^v are purely spatial (and symmetric) tensors. 

If := (t, x 1 ) are the spacetime coordinates adapted to the 
foliation of M introduced above, the line element is usually writ- 
ten in the so-called ADM form 

ds 2 := g flv dx fJ dx v = -a 2 dt 2 + y ij (dx i + 0dt)(dx j + fidt), (4) 

where the lapse function a and the shift vector /? (a purely 
spatial vector) are free gauge functions. In this adapted coor- 
dinate system the unit normal vector has components n M = 
(l/a,-P'/a) and = (-a, 0,). In the 3 + 1 formalism, the 
Einstein equations of Eq. (1) are split into a set of evolutionary 
equations for y t j and the extrinsic curvature Kjj 

d t7i j = -laKij + Dfij + Djpi, (5) 

d,K u = /3 k D k K u + K ik Dfi k + K jk Dfi k - DiD ja + 

a[R u + KKij - 2K ik K)} + 4na[y ij (S -E)-2S y ], (6) 

plus a set of constraints that must be satisfied at all times 

R + K 2 - KijK ij = \6nE, (7) 

Dj(K ij -Ky ij ) = SnS\ (8) 

named Hamiltonian and momentum constraints, respectively. 
Here D, := y^V^ is the covariant derivative with respect to the 3- 
metric y t j (so that D k yij = 0), Rij is the Ricci tensor, again with 
respect to y,j, R := R\ the corresponding Ricci scalar, K := K. is, 
the trace of the extrinsic curvature. As far as the fluid sources are 
concerned E := n„n Y T^, S 1 := -n^T^, and S ij := f^T^ 
(of trace S := Sj) are, respectively, the energy density, momen- 
tum density, and the stress-energy tensor as measured by the 
Eulerian observers. 

The constraints introduced above are notoriously difficult to 
be maintained in the numerical evolution of Eqs. (5) and (6), and 
two possible approaches can be followed. The most widely used 
one relies on hyperbolic formulations of the initial value prob- 
lem, the constraints are imposed only for the initial data and nu- 
merical errors are just monitored or damped during time evolu- 
tion (free evolution schemes). On the other hand, the constraints 
can be enforced at each timestep during the numerical simu- 
lation, leading to the so-called constrained evolution schemes, 



where the main idea is to maximize the number of elliptic equa- 
tions, usually more stable than hyperbolic equations. Moreover, 
in the steady state case the set of equations should easily re- 
duce to those used for stationary spacetimes and for the construc- 
tion of initial data. An example is the so-called fully constrained 
formalism (FCF) for asymptotically flat spacetimes in full GR 
(Bonazzola et al. 2004), which contains the widely used set of 
CFC elliptic equations as an approximation. In the following we 
will describe the general assumptions of conformal flatness and 
we will review the set of CFC equations. 

Let us start by applying a Lichnerowicz conformal decom- 
position 

7U := >!> : Jn. (9) 

assuming a flat background metric fj (time independent and not 
necessarily in Cartesian coordinates), where the conformal fac- 
tor satisfies iff = (y/f) 1/12 , with / := detfj. A second assump- 
tion is the condition of maximum slicing of foliations 

K = 0. (10) 

Under these assumptions, to be preserved during time evolution, 
the trace of Eq. (5) and its traceless part become, respectively 

d, In y 1/2 = 5,1mA 1 /6 = Dfi, (11) 

2aK u = DjfSj + Dfii - |(D^) r , 7 , (12) 

where d t y — if, and only if, Df? = and the extrinsic curvature 
can be expressed in terms of derivatives of the shift vector alone. 
The next step is to use covariant derivatives associated to the 
flat 3-metric fj, that will be indicated here by the usual nabla 
operator V, (and V k fj = 0). When Eq. (9) holds, it is possible to 
demonstrate that the Ricci scalar for jij (that for fj is zero) is 

R = -8i/r 5 Ai/f, (13) 

in which A := V,V' is the usual Laplacian of flat space. The 
above relation combined first with the Hamiltonian constraint in 
Eq. (7) and then with the trace of Eq. (6) provides the following 
two scalar Poisson-like equations for the conformal factor if/ and 
the lapse function a 

Ai// = - [2nE + iKijK'-i] if, 5 , (14) 

A(aiff) = [2n(E + 2S) + \KijK 1 '] aifr 5 , (15) 

where we still need to write Kij in terms of flat space derivatives 
of the shift vector. 

To the traceless extrinsic curvature A,y := Kij - \Kjij = K t j 
is then applied a conformal time-evolution rescaling 

K ij = i/T 4 A' 7 , 2oA' 7 := (Ly6)' 7 , (16) 

where the conformal Killing operator associated to the flat metric 
and applied to the vector /?' is defined as 

(Lj3) ij := V'fi + V j - l(V k j3 k )fj. (17) 

This scaling is also employed in the so-called conformal thin 
sandwich (CTS) approach to initial data. Moreover, it is quite 
natural, due to Eq. (12), and since within a conformally flat de- 
composition of the metric we have 

(LyfiyJ = ifr-\L/3r, (18) 

where L y indicates here the conformal Killing operator associ- 
ated to jij. On the other hand, in conformal flatness we also have 

DjK'j = if/- m Vj(if/ w K ij ), (19) 
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to be used with the above rescaling in the momentum constraint 
to find an equation for /?. 

Thanks to all the relations derived so far, the final set of CFC 
elliptic equations may be written in terms of the sources and of 
A'-* (containing a and first derivatives of /?) as 

Aif, = -[2nE + tfikfjt&jA") if/ 5 , (20) 

A(atff) = [2n(E + 2S) + IfikfjiAW] aif/ 5 , (21) 

A L /3' = I6naif/ 4 S 1 + 2if/ b A ij V jiaif/' 6 ), (22) 

where 

A L := Vj(L/3) iJ = A/3' + ±V'(V^), (23) 

is the so-called conformal vector Laplacian operator, associated 
to the flat 3-metric fj and applied to /3'. 

2.2. From CFC to XCFC 

A slightly different approach to the Einstein equations for 
asymptotically flat spacetimes has been recently presented 
(Cordero-Carrion et al. 2009). This involves a rewriting of the 
elliptical part of the FCF system for full GR through a differ- 
ent decomposition of the extrinsic curvature. Here we will just 
describe its conformal flatness approximation, leading to the so- 
called extended conformal flatness condition (XCFC) system of 
elliptic equations, improving on the CFC ones described in the 
previous section. The set of XCFC equations is our choice for 
the metric evolution in X-ECHO. 

The new approach still relies on the usual conformal de- 
composition in Eq. (9) and on the maximum slicing condition 
of Eq. (10), but the choice for the decomposition of the (trace- 
less) extrinsic curvature is different. We use here the momentum- 
constraint rescaling and the so-called York conformal transverse 
traceless (CTT) decomposition, first introduced for initial data, 
that is 

K ij = i//~ w A ij , A ij := (LW) ij + A^ T , (24) 

where the conformal Killing operator associated to the unknown 
vector W gives the longitudinal part of A'^ T , whereas Aj T 

is a transverse (V 7 Aj? r = 0), traceless (fjAj T = 0) tensor. 
Consistency between the CTS and CTT decompositions (notice 
that A'i = if/ 6 A'j) should require a non-vanishing A^ T . However 
it has been demonstrated that this quantity is even smaller than 
the non-conformal part of the spatial metric within the CFC ap- 
proach, and hence can be safely neglected on the level of the 
CFC approximation. Thus we set, as an additional hypothesis 

A^ T = => A ij = (LW) ij , (25) 

so that A'j is defined in terms of the auxiliary vector W alone. 
The latter is to be derived from the momentum constraint using 
Eq. (19), that is simply 

VjA ij = A L W' = 87ufr w S' , (26) 

to be added to the other CFC equations. 

The final augmented set of CFC elliptic equations, also 
known as XCFC equations, is then the following 

A L W* = 8nf ij Sj, (27) 

Aiff = -2nE if,- 1 - IfikfjiAW r\ (28) 

A(ar^) = [2n(E + 2S ) i/T 2 + IfafyA^ r 8 ] aif/, (29) 



A L = I6n ai//~ 6 f u Sj + 2A' 7 V/m/T 6 ), (30) 

where for convenience we have introduced rescaled fluid source 
terms of the form 

S j-.= if/ 6 S j, E:=if/ 6 E, S:=if/ 6 S, (31) 

and we remind that 

&H = y w J + V j w «' _ i(V k W k )fj. (32) 

Some comments and comparisons between the CFC and XCFC 
sets of equations are now due. 

- The unknown functions are now 8 rather than 5 (W, iff, a, 
/?), and this is reflected by the augmented number of ellip- 
tic equations (there is a new vector Poisson equation for the 
auxiliary variable W). 

- While in CFC all the equations were strongly coupled, here 
the equations can be solved hierarchically one by one, in the 
given order, since each right hand side just contains known 
functions or the variable itself (in the two scalar Poisson-like 
equations for if/ and aif/). 

- As we will see in the next sub-section, schemes for general 
relativistic hydrodynamics or MHD (like ECHO), given a 
metric in 3 + 1 form, actually evolve in time the conserva- 
tive variables y ll2 S j and y 1/2 E, rather than S' and E. Since 
if/ 6 = y l l 2 1 ' f x l 2 and f 1 ^ 2 is known and time-independent, 
the sources S j and E are basically known after each com- 
putational time-step without the need of an updated value 
of if/. This will be only needed to work out S = if/ 6 yijS lj , 
after the new value of if/ has been provided by Eq. (28) 
and the inversion of conservative to primitive variables has 
been achieved. Primitive variables are then updated self- 
consistently together with the new values for the metric, 
whereas this was not possible in CFC. In that case, one could 
either use Eq. (11) to derive a guess of the updated if/ (a 
method easily prone to both convergence problems and dis- 
cretization errors), or one is forced to iterate simultaneously 
over the metric solver (the whole CFC set) and the inversion 
routine for the primitive variables (typically itself a numeri- 
cal iterative Raphson-Newton method). 

- The last, and certainly not least, issue is related to the math- 
ematical nature of the scalar Poisson-like equations. In both 
cases we have a structure of the form 

Am = hu p , (33) 

where u is the generic variable (if/ or aij/), h is the generic 
source term, and p provides the exponent of the non-linearity 
(p — for a canonical Poisson equation). It can be demon- 
strated that the condition ph>0 implies that the solution u is 
locally unique. While this is always true in XCFC, since we 
have two contributions with p — - 1 and p = -7, both with 
h < 0, in Eq. (28), and one contribution with p — +1 and 
h > in Eq. (29), local uniqueness cannot be guaranteed for 
the CFC system, since Eq. (21) contains a term that certainly 
violates the requirement (the second one, due to the presence 
of a factor a^ 1 in A'- 7 ). 

2.3. The GRMHD equations and the ECHO scheme 

The equations for an ideal, magnetized, perfectly conducting 
plasma are 

V> M ") = 0, (34) 
= 0, (35) 
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V/*"' = 0, 



(36) 



respectively the continuity equation, the conservation law for 
momentum-energy, and the sourceless Maxwell equation. Here 
p is the mass density as measured in the frame comoving with 
the fluid 4- velocity u M , the total momentum-energy tensor is 



T» v = ph u^u v + pg>™ + F" A F vA - \{F M F AK )g^\ 



(37) 



with h = 1 + e + p/p the specific enthalpy, e the specific inter- 
nal energy, p = pip, s) the thermal pressure (provided by some 
form of equation of state, EoS). Moreover, F MV is the Faraday 
(antisymmetric) electromagnetic tensor, with the associated dual 
F*" v = \e^ vM F M , where e^ K = {-gY ll2 [pvAK\ is the space- 
time Levi-Civita pseudo-tensor = -(-g) l/2 []uvAK]), with 
g = detg^y and [[ivAk] is the alternating Levi-Civita symbol. 
The system is closed by Ohm's law for a perfectly conducting 
plasma, which becomes a constraint for a vanishing electric field 
in the frame comoving with the fluid 



F^Uy = 0. 



This basically replaces the Maxwell equation V^F^ 



(38) 
-J\ 

where the 4-current J v is a derived quantity as in classical MHD. 

In order to derive the GRMHD equations in 3 + 1 form, as 
employed in ECHO, we must decompose all 4-vectors and ten- 
sors into their spatial and temporal components on each slice 
of the time evolution. This can be easily achieved by using the 
unit normal vector introduced in the previous section 



puf := D(v" + n"), 
r" v := S" v + n^S v + S"n v + En^n v , 



F" v := n tl E v - E^n v 



B A n K 



n»B v - B»n v - e^Em 



(39) 

(40) 
(41) 
(42) 



where every new quantity is now purely spatial, as measured by 
the Eulerian observer. In particular, D := pY is the rest mass den- 
sity, v' is the fluid velocity, F := (1 - v,V)~ 1/2 is the usual Lorentz 
factor, whose definition is due to the condition u^u^ = -1. The 
stress-energy 3-tensor S' j , its trace S, the momentum density 
S', and the energy density E are the same quantities appearing 
in the 3 + 1 Einstein equations (this is why notations have been 
slightly modified with respect to the original ECHO paper), and 
are respectively given by 

S u = phY\vj + p Tl j - EiEj - BiBj + \(E k E k + B k B k )y u , (43) 
S = ph(Y 2 -l) + 3p + {{EiP + BiB'), 



S i =phY 2 v i + e ijk E j B k , 



(44) 
(45) 

E = phY 2 -p+ \{EiE* + B^). (46) 

The electric and magnetic fields as measured by the Eulerian 
observer are defined as E' := -n il y l v F iIV and B' := -n^F*^ , 
and due to the condition in Eq. (38), the electric field is a derived 
quantity precisely as in classical MHD 



E t = -e ijk v j B\ 



(47) 



where e ijk = y ' [ijk] (e' Jk = y ' [ijk]) is the Levi-Civita 
pseudo-tensor for the 3-metric y,- 7 -, and [ijk] is the alternating 
symbol taking values +1, -1, or 0. 

Thanks to the above decompositions, the GRMHD equations 
can be entirely rewritten in terms of purely spatial vectors, while 



retaining the original conservation form. We end up with a sub- 
set of fluid-like balance laws in divergence form 



3fU + diT' = S, 



(48) 



plus a magnetic sub-set with the induction equation in curl form 
and the associated divergence-free condition, to be preserved at 
all times during evolution 



d0 + [ijk]dj& k = 0, 30 = 0. 



(49) 



The set of conservative fluid variables and the set of associated 
fluxes are, respectively 



<U :=y 



1/2 



i ■- -,1/2 



r :=y 



(av'-p^D 
aS'j -0Sj 
aS'-ftE 



(50) 



whereas the set of source terms contain the derivative of the met- 
ric and thus the curvature effects 



1/2 





aS^djyik + SidjB' -Edja 



aK u S'J 



■S J d ja 



(51) 



Here Kij is the extrinsic curvature introduced in the previous sec- 
tion, whose evolution is directly provided by the Einstein equa- 
tions in the 3 + 1 formalism, together with y^, or may be given 
in terms of the derivatives of the metric terms. For a dynamical 
spacetime under the maximal slicing condition, from Eq. (12) 
we can write 

aK u S iJ = \pS ik djy ik +S\djP - |[r 1/2 5,<y 1/2 ^')]S, (52) 

and the same expression without the last term may be used for 
any stationary spacetime (Cowling approximation). As far as the 
induction equation is concerned, we have defined 



£':=y 1/2 B\ 



(53) 



+ e ijk B j B k = -y 1 ' 2 [ijk](av j -p)B k , (54) 

where the vector av' - B' is sometimes called transport velocity. 
The induction equation may be also written in the equivalent 
divergence form as 

d,& + 5,{y 1/2 [(av'' -B i )B j - B\av j -B j )]} = 0, (55) 

and the antisymmetric nature of the magnetic fluxes reflects that 
of the original electromagnetic tensor. Whatever is the adopted 
choice for the form of the induction equation, we have a final set 
of eight hyperbolic equations. Usually the corresponding vari- 
ables are named primitive variables, for example the set 



P:= [p,v\p,B l ] 



hT 



(56) 



for which <U = T l = T\P), and S = S(P), where for 

simplicity we consider here the augmented system with & in 

As discussed in the previous section, the inversion of the 
non-linear system It = ( U( e P) to recover the set of primitive 
variables is achieved through a numerical iterative scheme and 
requires the knowledge of the volume element y 1 ^ 2 , for consis- 
tency updated at the same time level as 1i. Moreover, as dis- 
cussed in Del Zanna et al. (2007), the conservative to primitive 
variables inversion is the most delicate part of a relativistic MHD 
code, as high Lorentz factor flows or strong magnetic fields (i.e. 
in a NS magnetosphere) may easily lead to errors in the values 
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of the conservative variables. The whole procedure can be re- 
duced to the solution of two coupled nonlinear equations, for any 
given EoS. In X-ECHO we leave complete freedom in the choice 
p = pip, s), however in the present paper, given the nature of the 
numerical tests proposed, we limit this choice to either an ideal 
y-law EoS 

pip,s) = iy-l)ps, (57) 
or to a polytropic law EoS 

pip) = Kp\ , sip) = -^—Kp-y-\ (58) 
J- 1 

where y and K are given constants. When the first option is used, 
the root-finding procedure solves for the variable x = v;v', ac- 
cording to the third method described in the ECHO paper, where 
the nested second equation (a cubic) for y = phT 2 can be either 
solved analytically or iteratively with an inner loop (and we typ- 
ically adopt this latter choice). Notice that when S ,B' = 0, or in 
general for a purely fluid simulation, the equation for y is linear 
and can be readily solved. On the other hand, for a polytropic 
law the energy equation becomes redundant, since all thermody- 
namical quantities are now functions of the density alone, and 
the overall inversion procedure simply reduces to a root-finding 
iterative solver for x. 

3. The XCFC solver for axisymmetric spacetimes 

The X-ECHO code for GRMHD in dynamical spacetimes is 
built upon the ECHO scheme, coupled to a novel solver for the 
XCFC equations described in Sect. 2.2. All the necessary defini- 
tions and equations have been provided in the previous section, 
here we specialize to the particular implementation of X-ECHO 
we are mostly interested in, that is under the assumption of ax- 
isymmetric GRMHD configurations, adopting spherical-like co- 
ordinates x' = (r, 6, ft) for the conformal flat metric. 

Thus, as a first step we assume in Eq. (9) the usual spherical 
coordinates 

/; v = diag(l,rVsin 2 0), (59) 

so that 

ds 2 = -a 2 dt 2 +ft 4 [idr+f3 r dt) 2 +r 2 id6+{?dt) 2 +r 2 sin 2 Oidft+p^dtf] 

(60) 

with f 1 / 2 = r 2 sin$ and y 1 ^ 2 = ft 6 r 2 sin9. Moreover, we are 
going to specialize to the axisymmetric case, thus the condition 

<9 = 0, (61) 

will be assumed throughout the paper. Within the flat metric 
fij, it is convenient to introduce the orthonormal basis e\ := 
ie f ,e g ,e^), with 

ef-=d r , e e :=r~ l d , e$ := (r sin0) _1 <9 , (62) 

for which fij = diag(l, 1,1), where a similar notation as in 
Bonazzola et al. (2004) has been assumed. Any generic vector 
X can be then expressed in the usual form as 

X:=X f e ? +X%+X%, (63) 

where the orthonormal (with respect to fij) components X' are, 
respectively 

X f := X r , X 6 := rX e , X^ := rsinflX*, (64) 



while the relation to covariant components still involves the 
function ft, since X[ = ft A fijXK As far as covariant derivatives 
are concerned, the change of basis allows one to use the V; op- 
erator of spherical coordinates. In particular, the Laplacian of a 
generic scalar function uir, 6) is 

Au = d 2 u + 2r~ l d r ii + r~ 2 id 2 e ii + coiOdgu), (65) 

whereas the orthonormal components of the conformal vector 
Laplacian are, respectively 

iA L X) f = AX f - 2r 2 iX r + d e X e + cot X s ) + ±<9 r (V • X), (66) 
iA L Xf = AX° + 2r 2 d g X f - (r sin ey 2 X" + \r~ l dgiV ■ X), (67) 

iA L Xf = AX* - (r sin 6»r V, (68) 
where the divergence of X is 

V • X = d r X f + r~ 1 (2X' r + d e X e + cot 9 x\ (69) 

precisely the formulae of vector calculus in spherical coordi- 
nates. 

The Poisson-like elliptic equations used in XCFC to com- 
pute the eight metric terms consist of two of scalar equations in 
the form of Eq. (33), for the variables u - if/ and u = aft, and 
two vector equations for the generic unknown vector X' = W 
and X' - P'. Due to non-linearity, the scalar equations are better 
solved iteratively for the quantity q n := u n — 1, which in both 
cases gives the deviation from asymptotic flatness ft — > 1 and 
a — > 1, where the new value at step n is computed using the pre- 
vious value at step n - 1 in the source term, until convergence 
is reached within some prescribed tolerance. Summarizing, the 
metric equations are expressed in one of the two generic forms 

Aq n = H n ^=hi\+q n - X ) p , (70) 

and, using the orthonormal basis introduced above, 

(A L *y = H 1 , (71) 

where h and H are generic scalar and vector source terms, to 
be provided by the XCFC equations, so that both the right hand 
sides above are known functions of (r, 6). 
, Numerical methods available for solving these elliptic par- 
tial differential equations (PDEs) can be divided into three main 
categories (see, for methods and discussions Grandclement et al. 
2001; Dimmelmeier et al. 2005; Grandclement & Novak 2009): 
direct inversion, full relaxation, spectral. Direct inversion codes 
are able to solve the complete system of CFC (or XCFC) at once 
using Newton-Raphson solver (or any other inversion technique) 
on the entire computational grid over which the metric solution 
is desired. They have very good convergence to machine accu- 
racy within few steps, but they suffer serious limitations: the ini- 
tial guess must be close enough to the solution to avoid con- 
vergence on local minima (instead of global minima); memory 
requirement for matrix allocation is typically very large (usually 
a sparse matrix is needed in the whole 2D domain); in general 
direct inversion schemes solve the metric on smaller grids that 
the one over which the fluid variables are evolved, and require 
interpolation between the two, with problems that may arise at 
boundaries. Full relaxation codes use SOR, multigrid or other 
relaxation techniques. The schemes are fast, they require little 
memory allocation, but usually suffer from poor convergence 
properties: the convergence is in general slow (compared to di- 
rect or spectral schemes) and might fail on the axis or at the cen- 
ter, due to the singular nature of some metric elements. Spectral 
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schemes decompose the set of CFC (or XCFC) equations using a 
combination of spherical harmonics (based on Legendre polyno- 
mials) in the angular directions and Chebyshev polynomials in 
the radial direction. This ensures a correct behavior on the axis 
and at the center even with a limited number of eigenfunctions, 
but they require specific grids and sometimes complex compacti- 
fications or multi-domain decomposition techniques with appro- 
priate boundary conditions for each multipole and each domain. 
The metric solver in X-ECHO uses a mixed technique. In the 
angular direction we decompose using spherical harmonics, to 
preserve the correct asymptotic form on axis. However, the set 
of ordinary differential equations (ODEs) obtained for each har- 
monic is then solved using direct inversion over the same radial 
grid used in the GRMHD code, with no need of interpolation or 
compactification. At second order accuracy in a finite difference 
discretization, the scalar equations reduce to the simple inversion 
of tridiagonal matrices (band diagonal matrices for the poloidal 
components of the vector Poisson equations), where appropriate 
solvers are fast, require little memory allocation, and typically 
converge with high accuracy. 

For the scalar Poisson-like Eq. (70) we then decompose, for 
each level n of iteration (here omitted), the unknown q as 



for Ci(r) alone. The new source terms are given by 



q(r,0) ■= [Ai(r)Yi(6)] , 



1=0 



where 



Y£9) = Y?(0) := J2£IPKcos0), 



(72) 



(73) 



with Pi the Legendre polynomial of degree I and the axisym- 
metry condition has been imposed (m = 0). The PDE, where 
the Laplacian is provided in Eq. (65), may be then split into the 
series of radial ODEs for each harmonic I 



d 2 A, 2 dA t 
dr 2 r dr 
where the new source term is 



H,(r) := j)H(r, 



0)Y,(6)dQ., 



(74) 

(75) 
to 



with dQ. = 2n sin 6d6 and the integral running from 6 
6 - n due to axisymmetry. 

As far as the vector Poisson equation in Eq. (71) is con- 
cerned, the unknown vector X is first decomposed into vector 
spherical harmonics, that is, in the axisymmetric case 

oo 

X(r, 6) := £ [A,(r)Y,(e)e f + B ^Yj (6)e d + QW^] , 

/=o 

(76) 

where Y\ := dY\jdQ. As it is apparent from the operators in 
Eqs. (66-69), the set of equations split into a series of ODEs 
with, for each harmonic /, a coupled poloidal part for the radial 
functions A/(r) and Bi(r) 



4 d 2 A, 8 / dAi \ 1(1+1) ( A rdB, 1 



d 2 B, 2dB, 41(1+1) n 1 dA t 



+ r dr 



3r 2 



3r dr 3r z 



dr 2 

and a toroidal part 

d 2 C, 2 dC, 1(1+1) 



dr 2 r dr 



— --^-Q=H 



Hi 
(77) 

(78) 



(79) 



H[(r) :- 



H r (r, 6)Yi(6)d£l, 



H (r) :-- 



1 



1(1+1) 



flf(r):= 



1 



/(/+ 1) 



§w 



(r,9)Y' l (9)dn, 



(r,6)Y',(6)d£l. 



(80) 



(81) 



(82) 



These integrals, as well as that in Eq. (75), are computed in 
X-ECHO by using Gaussian quadrature points, so the original 
source terms must be first interpolated on the required locations. 

As anticipated, we use finite differences and a second order 
approximation for first and second spatial derivatives, so that, for 
each harmonic /, the two coupled poloidal equations are reduced 
to the inversion of a sparse matrix of bandwidth 7, whereas the 
toroidal equation leads to a tridiagonal matrix (standard open 
source routines are employed). Boundary conditions at r — 
and r = r max are given by the parity and asymptotic properties of 
the multipole corresponding to the harmonic /. In particular, we 
assume that at the center A/(r) has parity (-1)', whereas Bi(r) and 
Ci(r) have parity (-1)' +1 , and at large distances from the central 
sources the multipoles are forced to decay as r~' (/+1) . 

4. Axisymmetric GRMHD equilibria: the XNS code 

One of the most common application of the fully constrained 
formalism is the search of self-consistent (axisymmetric) equi- 
librium configurations (i.e. fluid quantities and metric) for com- 
pact relativistic stars (we will simply refer to these objects with 
the term NS) a typical case of initial data problem in GR (Cook 
2000; Stergioulas 2003; Gourgoulhon2010). Several codes have 
been presented in the years addressing this issue (Komatsu 
et al. 1989a,b; Cook et al. 1994; Stergioulas & Friedman 1995; 
Nozawa et al. 1998; Bonazzola et al. 1998; Kiuchi & Yoshida 
2008), and all of them, despite the different approaches and up- 
grades (e.g. differential rotation, toroidal magnetic field), adopt 
the so-called quasi-isotropic coordinates. Under the conditions 



8, = 0, d, = 0, 



(83) 



we write here the corresponding line element in the form 

ds 2 = -a 2 dt 2 + tfr 4 (dr 2 + r 2 d6 2 ) + R 2 (dcf> - codt) 2 , (84) 

which resembles that of the CFC metric for axisymmetric space- 
times of Eq. (60) in spherical coordinates and reduces to it when 



R = limine, 



(85) 



where only in this case the function iff in Eq. (84) recover the 
meaning of conformal factor. In general, we can think the metric 

1 II 

function R :- to be a sort of generalized cylindrical radius, 
whereas a> := is the intrinsic angular velocity about the sym- 
metry axis of the zero angular momentum observers (ZAMOs: 
Bardeen et al. 1972), that are the normal (Eulerian) observers 
for an axisymmetric spacetime. When oj = the spacetime is 
spherically symmetric and the two metrics both reduce to that in 
isotropic Schwarzschild coordinates. In the following we briefly 
re-derive the GRMHD equilibrium condition (a Bernoulli-like 
integral) for purely toroidal flows and magnetic fields in quasi- 
isotropic coordinates, more general derivations can be found in 
(Kiuchi & Yoshida 2008) or (Komissarov 2006), this latter work 
applied to magnetized tori around rotating BHs (see also the final 
numerical test in Del Zanna et al. 2007). 
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The only non vanishing components of the spatial (Eulerian) 
3-vectors v' and B' are the azimuthal ones, and the corresponding 
moduli are 

v : = O v^) 1/2 = Rv* = cT l R{Cl - co), (86) 

B := (B B ) 1/2 = RB^, (87) 

where we have used the 3+1 relations = F(v^ + oj/a), 
u' = F/a, and we have defined CI := u^/u' = d<p/dt, the an- 
gular velocity of the fluid seen by an observer at rest at infinity 
(a function of r and 6 for a differentially rotating NS). The equi- 
librium condition we are looking for can be directly derived from 
the 3 + 1 conservative form of the GRMHD equations, in which 
the only non vanishing contribution is due to the two poloidal 
components of the momentum conservation in Eq. (48) 

y- 1/2 dj[ 7 ' /2 a(p + \B 2 )] - a(p + \B 2 )\rdf/u 

- a(phY 2 v 2 - B 2 )\R- 2 djR 2 + phT 2 RvdjCo 

+ (phY 2 -p + \B 2 )dja = 0, (88) 

where j is either r or 9 and we remind that the electric field van- 
ishes. Recalling that \y ll djyu = y~ 1/2 djy 1/2 and differentiating 
the definition of v, after a few algebraic steps the following con- 
dition can be found: 

dip Y 2 v 2 d,Cl dj(a 2 R 2 B 2 ) 

-!f + djlna - djlnT + - 1 + -£——L = 0, (89) 

ph 1 1 CI — a) 2a 2 R 2 ph 

Integrability of this equation demands the following conditions: 

- A barotropic EoS, p = p(ph). The simplest choice and the 
most common assumption is a polytropic law 



p = Kp l+l/n => ft = 1 + (n + l)Kp l/n , 



(90) 



where n is the polytropic index (the corresponding adiabatic 
index is y = 1 + 1/n). 
- The quantity F := u^u' = Y 2 v 2 /(Cl-co), related to the specific 
angular momentum t := -u^/u,, is a function of the angular 
velocity CI alone. A commonly adopted differential rotation 
law (e.g. Stergioulas 2003) is 



, R 2 (C1 - co) 

F(Cl) = A 2 (Q C - Q) 



a 2 -R 2 (Cl-co) 2, 



(91) 



where Cl c is the central angular velocity and A is a measure 
of the differential rotation rate. For uniform rotators Cl = Cl c 
(and A — > oo), and this contribution can be excluded from 
the Bernoulli integral. 
- A sort of magnetic barotropic law, where aRB is a function 
of a 2 R 2 ph. The simplest choice is a magnetic polytropic law 



aRB = K m (a 2 R 2 ph) m , 



(92) 



where m is the magnetic polytropic index, with m > 1 
(Kiuchi & Yoshida 2008). 

Using the above prescriptions we can easily derive the final 
GRMHD Bernoulli integral 

In L + i n ± _ i„ r - -{Cl c - Q) 2 + J!^L( a 2 R 2 ph) 2m - 1 =0, 
h c a c 2 2m- I 

(93) 

where again with the subscript c we indicate values at the center. 

The above derivation is exact for quasi-isotropic coordi- 
nates and obviously applies also to the CFC sub-case. In the 



non-magnetized case, it has been demonstrated that, even for 
rapid rotators close to the mass shedding limit (see Sect. 5.3), 
solutions obtained with the RNS numerical code (Stergioulas 
& Friedman 1995) show very little deviations from a CFC 
metric, so in principle one would expect the CFC limit to 
provide a reasonably good approximation of the correct solution 
for compact relativistic stars. Given the above mentioned 
computational advantages in the solution of the XCFC system 
of equations with respect to the original set of coupled PDE of 
fully constrained schemes, our choice has been to build a novel 
numerical code (written in Fortran9Q), which we name XNS 
and that can be freely downloaded at 

http : // sites . google . com/site/niccolobucciantini/xns 

This takes advantage of the same XCFC metric solvers 
developed for the X-ECHO scheme described in the present 
section. A comparison between equilibria found with XNS and 
RNS is presented in Sect. 5.3, together with a discussion of the 
results. 

XNS employs a self-consistent method in the search for the 
axisymmetric equilibrium solutions of relativistic compact stars, 
in the presence of differential rotation and a purely toroidal mag- 
netic field, thus metric terms and fluid-like quantities are derived 
at the same time. Given the values of the six free parameters 
K, n, A, Cl c , K m , m, plus a guess for the central density p c , the fol- 
lowing steps are taken: 

- A starting guess for the CFC metric terms (a, if/, co) is pro- 
vided from the previous step, with R given by Eq. (85). The 
first time, the spherically symmetric Tolman-Oppenheimer- 
Volkoff (TOV) solution in isotropic coordinates (Tolman 
1934), for the metric of a non-rotating and non-magnetized 
NS with a central density value p c , is computed through a 
shooting method for ODEs. 

- On top of these metric terms and for each grid point, CI is de- 
rived by inverting Eq. (91), then v (and T) can be determined 
from Eq. (86). Finally, h c is a known function of p c and the 
Bernoulli integral in Eq. (93) is solved via a Newton method 
to find the local values of h and p, allowing us to determine 
also the magnetic field strength from Eq. (92). 

- Combining the updated fluid quantities with the old metric, 
the new conserved quantities E, § j are derived. 

- The set of XCFC equations is solved for this set of con- 
served variables, and a new metric is computed. Here only 
the azimuthal components for the vector Poisson equations 
are treated, in order to find and then B^. We remind here 
that the conservative to primitive variables inversion must be 
enforced between the solutions to the two scalar Poisson-like 
equations, namely after the new value of if/ is found and be- 
fore that of a. 

These steps are repeated until convergence to a desired tolerance 
is achieved. However, a word of caution is in order here. Since 
XNS is based on the XCFC metric solver, that works on the con- 
servative variables densitized with if/ 6 , convergence is actually 
enforced on the central value of the quantity D := if/ 6 D = iff 6 p c . 
Therefore, given that the final conformal factor tf/ for the self- 
consistent 2D equilibrium may be quite different from that de- 
rived from the radial TOV solution at the first step, the final value 
of p at the center is expected to differ from the parameter p c in the 
Bernoulli integral. If one wants to find an equilibrium converg- 
ing exactly to a central density p c , an additional overall iterative 
loop is needed. 
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Before concluding the section, some remarks are in order. 
One might question the choice of a purely toroidal field ver- 
sus a more realistic configuration including a poloidal compo- 
nent. However, equilibrium with poloidal fields is only possible 
for uniform rotators, with magnetic field fully confined within 
the star. It is well known (Goldreich & Julian 1969) that any 
poloidal magnetic field extending outside a rotating NS will lead 
eventually to a spin-down of the same, even for a dipole aligned 
with the rotation axis, unless the star is surrounded by a true 
vacuum. A small charge density of order of 10 20 cirT 3 pairs is 
enough to break this assumption, even in the case of millisec- 
ond rotators with B ~ 10 17 G, and the timescale to replenish 
an evacuated magnetosphere is of order of the rotation period. 
For such strong magnetic fields, affecting the overall equilib- 
rium of the NS and providing a non negligible contribution to 
the global stress-energy tensor in the Einstein equations (lower 
fields have little dynamical effects and can be easily treated as 
perturbations) the problem becomes particularly severe because 
the typical spin-down time for millisecond rotators can be as 
short as 50 ms. Weaker magnetic field can lead to much longer 
spin-down times, and those configurations might be considered 
<7MG«/-equilibrium cases, at least for the time of a typical numer- 
ical run (a few thousands M). However they are of little interest, 
as stated above. 

One might also question if purely toroidal configurations 
are stable, and, if not, what is the growth rate of instabilities. 
A recent study of the stability properties of neutron star with 
strong toroidal field has been presented by Kiuchi et al. (2008). 
Their results show that non-rotating neutron stars with a mag- 
netic polytropic law with m > 2 are always unstable, and that 
(uniformly) rotating systems are stable only if their kinetic en- 
ergy exceed by at least a factor 5 the magnetic energy. However, 
even if their conclusions about stability are supported by numer- 
ical simulations, it must be pointed out that their criterion is de- 
rived analytically only in the Newtonian regime, where magnetic 
field and spacetime metric are not coupled. Moreover, their full 
GR results do not consider intermediate cases with 1 < m < 2, 
thus further investigation is certainly needed. 

5. Numerical results 

We present here a set of numerical tests of the X-ECHO scheme. 
Standard HD/MHD tests in a static background metric (Cowling 
approximation) have been already presented elsewhere both for 
a flat metric (Del Zanna & Bucciantini 2002; Del Zanna et al. 
2003) and for a given stationary curved spacetime (Del Zanna 
et al. 2007). As discussed in the introduction, the performances 
of the ECHO code in astrophysical scenarios involving a variety 
of different conditions like strong shocks, relativistic outflows, 
strong gravity, and highly magnetized systems, have also been 
already assessed in numerous papers. For these reasons, the re- 
sults presented here focus mostly on evaluating the quality of 
our novel metric solver, and the performances of the HD/MHD 
ECHO algorithms when coupled with a dynamical spacetime. 
Moreover, the original ECHO scheme was designed for high- 
order accuracy, whereas the metric solver implemented in X- 
ECHO is formally only second order in space. Given that the 
performances of high-order reconstruction techniques has al- 
ready been tested (Del Zanna et al. 2007), for simplicity, we 
have decided to limit our set-up to second order accuracy, both 
in space and time, which is always a good compromise between 
efficiency, accuracy, and robustness. 

As discussed in the introduction, only in recent years sta- 
ble numerical schemes for GRMHD in dynamical spacetimes 



have started to appear. This, together with the intrinsic degrees 
of freedom in the choice of gauge and coordinate systems, typi- 
cal of GR, have resulted in a lack of well defined, agreed upon, 
standard numerical tests (in the spirit of what shock-tube prob- 
lems are in flat spacetime or accretion/outflows solutions are for 
a stationary curved metric). While in ID a few problems have 
emerged as standard benchmarks, this cannot be said about mul- 
tidimensional cases yet, also because of a lack of analytical so- 
lutions for fully multidimensional dynamical problems. Some of 
the tests that we have selected here, have been already published 
in the literature, using different numerical schemes, both fully 
constrained (Dimmelmeier et al. 2002; Stergioulas et al. 2004; 
Dimmelmeier et al. 2006; Cordero-Carrion et al. 2009) and hy- 
perbolic (Font et al. 2002; Bernuzzi & Hilditch 2010), so that 
we can validate our code against existing results. Some other 
have been done in the perturbative regime and compared with 
the results of linear theory, and when possible also with previous 
numerical simulations. However, in an attempt to evaluate the 
performances under different conditions, we also present some 
novel cases. 

Unless otherwise stated, in all our numerical tests we will 
use a Courant number of 0.4, an ideal gas EoS with y = 2 (cor- 
responding to a polytropic index n — 1 for the initial data in 
XNS), and we will solve the full GRMHD system, including 
the equation for the total energy density E, often neglected in 
isentropic tests for a given polytropic law. For the approximate 
Riemann solver, we use here for the first time the HLLC solver 
by (Mignone & Bodo 2006), never applied before to GRMHD 
studies in curved, evolving spacetimes, to our knowledge. This 
choice is imposed by the sharp transition between the rotating 
NS and the external atmosphere, since the two-wave solver HLL 
is found to be too diffusive on the contact discontinuity. 

Given that the astrophysical problems we are mostly inter- 
ested in, and toward which the X-ECHO scheme has been devel- 
oped, involve the ability of simultaneously handle the high den- 
sity NS and any low density outflow/atmosphere/magnetosphere 
that might surround it, in almost all of our tests we have in- 
cluded in the domain an extended atmosphere, which is left free 
to evolve and respond to the evolution of the NS. We are aware 
that in the literature a common practice is to reset floor values 
outside the NS, but we believe that this procedure may in prin- 
ciple lead to violations of conservation properties of the scheme 
at the NS surface. For the same reasons, as stated above, simula- 
tions have been performed using a ideal gas EoS, which allows 
us to handle systems where matter in different thermal condi- 
tions (a cold NS versus a hot atmosphere) is present. 

Reconstruction at cell boundaries is achieved for simplicity 
through a monotonized-central (MC) algorithm, though the other 
choices described in the appendix of (Del Zanna et al. 2007) are 
also possible. The XCFC metric solver is invoked every 10 steps 
of the HD/MHD Runge-Kutta evolution scheme. In 2D runs we 
use 50 Gaussian quadrature points and 20 spherical harmonics. 
With these settings, we have managed to make comparable the 
computational times taken by the metric solver, usually slower, 
and by a single Runge-Kutta cycle of the fluid solver. Thus, solv- 
ing XCFC at every fluid step will only double the overall times, 
but no significant improvement has been noticed in the results. 
Finally, grid spacing will be constant both along the radial direc- 
tion r and the polar angle 9, so the number of points is enough to 
specify the grid in each direction. 
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Fig. 1. Evolution of a stable TOV solution in spherical symmetry and 
isotropic coordinates. The upper panel shows a comparison between 
density in the initial solution (solid line) and the result at f max = 1500 
(diamonds). For clarity the result a f max is shown every 5 points. The in- 
sert shows the residuals. The spike at r ss 8 is due to diffusive relaxation 
at the NS boundary. The middle panel shows the relative variations in 
time of the central density. The lower panel shows the Fourier transform 
of the the central density. Solid line and diamonds indicate the power 
of the Fourier series in arbitrary units. The vertical markers indicate the 
frequency of known eigenmodes. The frequency resolution of our time 
series is ~ 150 Hz. 



5.1. Stability of a TOV stable radial solution 

Our first test consists in the evolution of a stable ID radial 
NS configuration. We adopt, as initial condition, a solution of 
the TOV equations in isotropic coordinates, corresponding to 
a polytropic gas with K = 100, n = 1, and central density 
p c = 1.280 x 10~ 3 , a model also known as AO (AU0) or B0 
(BU0) in the literature (Font et al. 2000; Stergioulas et al. 2004; 
Dimmelmeier et al. 2006). This corresponds to a NS extending 
to a radius r = 8.13. It is possible to show that this star lies 
on the stable part of the mass-radius curve, so we expect the 
code to be able to maintain this configurations for times longer 
that their typical sound crossing time (~ 0.5 ms). Outside the 
star we assume, at the beginning of the run, a low density and 



hot (p ^ 10 -7 , pip ^ 0.2) atmosphere in hydrostatic equilib- 
rium {ah = const), in pressure balance at the surface of the NS. 
Contrary to previous treatments, where the atmosphere was reset 
at every time step to keep it stationary, we leave it free to evolve 
(collapse or expand) in response to the NS oscillations. Given 
its low density, the atmosphere has negligible feedback on the 
star. The simulation is performed using 625 grid points in the 
radial domain r = [0, 20], corresponding to a star resolved over 
250 points. The evolution is followed for a time f max = 1500 
corresponding in physical units to ^ 7.5 ms. 

Fig. 1 shows a comparison between the initial density pro- 
file and that at f max , together with a plot of the central density 
p c as a function of time. Relative variations of density in the NS 
interior are of order of 10~ 4 , with major deviations only at the 
contact discontinuity of the NS surface, due to diffusion over the 
much lower density atmosphere. This triggers the natural vibra- 
tion modes of the NS, that are the observed fluctuations, which 
are the natural outcome for this physical system. Notice that the 
central density, plotted in the insert of Fig. 1, shows fluctuations 
of order of 10~ 3 at most, but no sign of any secular trend. This is 
due to the large number of points over which the star is resolved. 
The slow damping of the oscillations, from 10~ 3 to a few 10~ 4 , 
is due to the thermal dissipation associated to the use of an ideal 
gas EoS and to the numerical viscosity of the scheme, whereas 
it is not present if a polytropic EoS is used (see next Sect. 5.2 
for a comparison between the two EoS). At the beginning of the 
simulation we observe a relaxation of the central density to a 
value which is ~ 2 x 10~ 4 lower than the initial condition, prob- 
ably because of discretization errors. However the average value 
seems to remain constant at later times. This shows the ability 
of the code to maintain a stable equilibrium, even for several 
(~ 10) sound crossing times. It also shows that the presence of 
a dynamical atmosphere causes no problem for the stability of 
the TOV solution. On the contrary, the atmosphere itself seems 
to be quite stable, as shown by the fact that its density changes 
by a factor smaller than 1%. 

In the bottom panel of Fig. 1 we plot a Fourier transform 
of the central density in time. Markers indicate the positions of 
the known eigenmodes. This is a test of the performance of the 
code in handling a dynamical spacetime, at least in the linear 
regime, for small perturbations. The values of the eigenmodes, 
the fundamental in particular, are quite different in the Cowling 
approximation where the metric is kept fixed in time (Font et al. 
2002). It is interesting to note also that no initial perturbation has 
been introduced, and that the oscillations of the star are only due 
to the relaxation of the initial conditions and possible round off- 
errors. Indeed, the large power present in the higher frequency 
modes suggests that the initial excitation is confined to small 
scales, most likely at the surface of the star. The presence of a 
freely evolving atmosphere does not seem to affect the frequency 
of the modes, at least within the accuracy of our temporal series. 

5.2. Migration of an unstable TOV radial solution 

A genuinely non-trivial situation in the fully non-linear regime, 
involving large variations of the metric and fluid structure, is the 
migration of an unstable ID TOV solution. Following (Font et al. 
2002; Cordero-Carrion et al. 2009; Bernuzzi & Hilditch 2010) 
we select a solution of the TOV equations in isotropic coordi- 
nates, corresponding to a polytropic gas with again K = 100, n = 
1, and central density p c = 7.993 x 10~ 3 . This corresponds to a 
star extending to a radius r = 4.26, on the unstable part of the 
mass-radius curve. The external atmosphere is set as in the pre- 
vious test. Due to truncation errors and the initial relaxation of 
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Fig. 2. Evolution of the central density for the migrating unstable TOV 
solution in spherical symmetry and isotropic coordinates: the solid line 
is the evolution for an ideal gas EoS, the dotted line is the evolution for 
a polytropic EoS. The horizontal dashed line indicates the density of the 
stable TOV solution corresponding to the same mass p c = 1.650 x 10~ 3 . 



the NS/atmosphere transition, this configuration migrates to the 
stable branch. This evolution causes large amplitude pulsations, 
and the formation of a shock between the outer mantle and inner 
core of the star, where part of the kinetic energy is dissipated 
into heat. During the evolution the star expands to quite large 
radii. The run is done using 900 grid points in the radial domain 
r = [0,30]. A small fraction of the stellar mass (smaller than 
~ 10~ 3 ) is lost at the outer radius during the first bounce. The 
evolution is followed for a time f max = 1500 corresponding in 
physical units to ^ 7.5 ms. 

Fig. 2 shows the evolution of the central density in time. 
Results agree with what has been previously presented, both 
concerning the amplitude of the fluctuations, the value of the 
density at the first minimum and at the first and second maxima, 
the frequency of the oscillations, their non-sinusoidal shape, and 
the asymptotic value of the central density with respect to the 
expected value for a stable configuration with the same mass. As 
already noted (Font et al. 2002), the lower average central den- 
sity at later time (with respect to that of a stable model with the 
same mass), is due to shock heating during the bounce phase, 
that changes the thermal content of the star. We have repeated 
the same simulation using a polytropic EoS for the NS (while the 
ideal gas EoS is still used for the free atmosphere). Results are 
shown in Fig. 2 (dotted line). Is is evident that for a polytropic 
EoS the oscillations are not damped, and the system appears to 
converge to the correct asymptotic value of the central density. 
Comparison with previous results (Font et al. 2002; Cordero- 
Carrion et al. 2009) agrees both in the amplitude of the oscil- 
lations and in their phase shift with respect to the case with an 
ideal gas EoS. 

5.3. Accuracy of uniformly rotating 2D equilibria 

Following Sect. 4, we have built a series of 2D equilibrium mod- 
els for rigidly rotating neutron star in the XCFC approxima- 
tion for the metric. The computational grid covers the domain 
r = [0,20] and 8 = [0, n], with 650 zones in radius and 100 
zones in angle. In this section we compare them with equivalent 
models built using the publicly available code RNS (Stergioulas 
& Friedman 1995; Nozawa et al. 1998; Stergioulas 2003), which 
for us constitutes a reference benchmark. RNS is an accurate 
solver for non-magnetized, uniformly rotating NS configurations 



Table 1. Comparison between RNS and XNS for rigidly rotating com- 
pact stars. Gravitational masses are in units of Mq, r e and r p are the 
equatorial and polar radii. 
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RNS 


XNS 


RNS 


XNS 


RNS 


XNS 


BU0 


0.000 


1.400 


1.400 


8.13 


8.13 


1.00 


1.00 


BUI 


1.075 


1.432 


1.433 


8.33 


8.34 


0.95 


0.95 


BU2 


1.509 


1.466 


1.468 


8.58 


8.56 


0.90 


0.90 


BU3 


1.829 


1.503 


1.505 


8.82 


8.85 


0.85 


0.85 


BU4 


2.084 


1.543 


1.455 


9.13 


9.15 


0.80 


0.80 


BU5 


2.290 


1.585 


1.586 


9.50 


9.52 


0.75 


0.75 


BU6 


2.452 


1.627 


1.629 


9.95 


9.98 


0.70 


0.70 


BU7 


2.569 


1.666 


1.667 


10.51 


10.53 


0.65 


0.65 


BU8 


2.633 


1.692 


1.693 


11.26 


11.30 


0.60 


0.60 


BU9 


2.642 


1.695 


1.698 


11.63 


11.60 


0.58 


0.58 



in quasi-isotropic coordinates, for which we recall that in 2D 
R 2 + if/ 4 r 2 sin 2 8 in Eq. (84). This is both a test of the quality of 
the XCFC approximation, for genuinely non-spherically sym- 
metric systems, and an evaluation of the performances of our 
metric solver. We have selected the BU series of models from 
the literature (Stergioulas et al. 2004; Dimmelmeier et al. 2006), 
with a fixed central density p c = 1.280 x 10~ 3 , the usual poly- 
tropic law with K = 100, n = 1, but different values of the uni- 
form rotation rates. Tab. 1 presents the models and compares 
some global properties derived using RNS with those derived 
using our implementation of the XCFC solver in XNS. Results 
agree within the accuracy of the models themselves (Ar = 0.03). 
The model BU9 represents an equilibrium at the mass shedding 
limit, and it is particularly sensitive to accuracy and round-off 
errors. Indeed, when solving the Bernoulli condition Eq. (93) to 
derive the matter distribution, in the case of model BU9, we had 
to impose the further condition p = for r > 1 1 .6, to avoid 
unbounded solutions. 

In Fig. 3 we show the comparison between the model BU8 
(we already pointed out the problems for model BU9) derived 
using our solver for the CFC metric and the solution of RNS. It is 
clear that the CFC approximation provides a good description of 
the matter and fluid properties of the equilibrium configuration, 
throughout the entire star. The densities along the polar axis and 
at the equator, together with the velocity profile v := (v^v*) -1 ^ 2 , 
are all well reproduced. In terms of the metric coefficient, we 
see that the discrepancy is of order 10~ 3 , and is comparable with 
the level of non conformally flatness, defined as the difference 
between iff and [R/(r sin 8)] 1 ^ 2 in Eq. (84). Somewhat larger de- 
viations, of order a few 10~ 3 in the NS interior, are characteris- 
tic of the shift ffi ', increasing at larger radii where the value of 
the shift approaches zero, analogously to what was observed by 
Dimmelmeier et al. (2002). 

5.4. Stability of uniformly rotating 2D equilibria 

In the same spirit as for the numerical test presented in Sect. 5.1, 
we show here the result of a time evolution of two equilibrium 
configurations, BU2 and BU8. They correspond to a mildly ro- 
tating star and to a case of rapid rotation, respectively. The ini- 
tial condition are derived according to the recipe given in Sec. 4 
in conformally flat metric, and their accuracy has already been 
discussed in the previous section. The domain, r — [0, 16] and 
8 = [0,n], is spanned by 200 zones in the radial direction 
and 100 zones in angle. The star is resolved on average over 
about 100 radial zones. The evolution is followed for a time 
fmax = 1500, corresponding in physical units to ^ 7.5 ms. As 
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Fig. 3. Comparison between RNS and XNS solutions, for model BU8. 
The upper panel shows the fluid quantities. The solid (dotted) line rep- 
resents the density at the equator (polar axis) derived using RNS, both 
normalized against p c . The dashed line is the profile of the rotational 
velocity module v, normalized to its maximum (0.37462). Diamonds, 
crosses and pluses represent values of the same quantities as derived us- 
ing XNS, where for clarity we report one symbol every 5 radial points. 
The lower panel shows the residual of various metric terms at the equa- 
tor. The green solid line is the relative error between the lapse a com- 
puted with XNS and RNS. Dashed magenta and dotted blue lines rep- 
resent the relative error between the conformal factor if/ of the CFC 
metric with respect to that in quasi-isotropic coordinates and the quan- 
tity [R/(rs'md)] >12 . The dot-dashed red line represents the difference 
between if/ and [R/(rsin0)] 1/2 both computed in quasi-isotropic coordi- 
nates, and can be considered a measure of the non conformal flatness of 
the RNS solution. 



in the previous tests, we initialize, outside the star, a hot and low 
density atmosphere, initially in equilibrium and then allowed to 
freely evolve in time. As was already noted in the previous ID 
cases, even in 2D the presence of this freely evolving atmosphere 
has negligible dynamical effects on the star, due to its low den- 
sity, and there is no need of enforcing a reset to the initial values. 

For the model BU2, Fig. 4 compares the equatorial and axial 
densities together with the module of the rotational velocity v, 
at time t — and t = f max . The evolution of the central density 
is also shown. The oscillations of order or 10~ 3 are due to the 
excitation of the stellar eigenmodes by initial round-off errors 
and relaxation at the boundary between the star and the atmo- 
sphere. No initial perturbation was introduced into the model. 
Changes in the value of the equatorial density at the end of the 
run, with respect to the initial values, are of the same order of 
the amplitude of the oscillations that are excited by the initial 
relaxation. A small secular drift of order ~ 5 x 10~ 4 in the av- 
erage value of the central density is visible by the end of the 
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Fig. 4. Evolution of a stable BU2 solution. The upper panel shows a 
comparison between the initial values (solid lines) of equatorial (green) 
and axial (blue) densities, and equatorial rotational velocity v (ma- 
genta), against the value of the same quantities at at f raax = 1500 (dashed 
lines). The densities are normalized to the central initial value p c , the 
velocity to its maximum initial value (0. 1619). The insert shows the rel- 
ative difference between the equatorial densities as a function of radius. 
The spike at r ~ 8 is due to diffusion at the NS surface, also partly 
visible in the velocity profile. The middle panel shows the variation in 
time of the central density. Note that there is a hint of a secular trend of 
increasing density, with an average increase during the whole run of a 
few 10~ 4 . The bottom panel shows the Fourier transform of the the den- 
sity (green solid line), v r (red dotted line), v* (magenta dashed line), and 
v e (blue dot-dashed line), at the point r = 3.0, 6 = 45° of model BU2. 
The vertical markers indicate the frequency of known eigenmodes. The 
frequency resolution of our time series is ~ 150 Hz. 



run. It is due to the lower number of points in the radial direc- 
tion over which the star is resolved (~ 100), compared to the ID 
case (~ 250) where no drift was observed. In the bottom panel 
of Fig. 4, we plot a Fourier transform of various fluid quantities 
in time. The quantities are all measured at a selected location 
inside the star: r = 3.0 and 8 = n/4. Normal mode analysis is 
beyond the scope of this paper, so here we just present a sim- 
ple single-point analysis. Selecting other points in the star does 
not change the location of the peaks in the power series (though 
it affects their relative amplitude). Given that no perturbation is 
introduced in the initial data, modes are differently excited, de- 
pending on the initial relaxation (certain modes can in principle 
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Fig. 5. Evolution of a stable BU8 solution. The upper panel shows a 
comparison between the initial values (solid lines) of equatorial (green) 
and axial (blue) densities, and equatorial rotational velocity v (ma- 
genta), against the value of the same quantities at at f raax = 1500 (dashed 
lines). The densities are normalized to the central initial value p c , the 
velocity to its maximum initial value (0.3499). The insert shows the rel- 
ative difference between the equatorial densities as a function of radius. 
The spike at r ^ 10 is due to diffusion at the NS surface, also partly 
visible in the velocity profile. The middle panel shows the variation in 
time of the central density. Note that there is a hint of a secular trend of 
increasing density, with an average increase during the whole run of a 
few 10~ 3 . The bottom panel shows the Fourier transform of the the den- 
sity (green solid line), v'~ (red dotted line), v* (magenta dashed line), and 
v" (blue dot-dashed line), at the point r = 3.0, 6 = 45° of model BU8. 
The vertical markers indicate the frequency of known eigenmodes. The 
frequency resolution of our time series is ~ 150 Hz. 



even not be excited). Markers indicate the positions of the known 
I — 0, I — 2 and I = 4 eigenmodes (Dimmelmeier et al. 2006). 
This is a test of the performance of the code in handling an ax- 
isymmetric dynamical spacetime, at least in the linear regime for 
small perturbations. The fundamental F and first overtone Hi of 
the I = mode are correctly recovered, together with the 2 p\ for 
the dipolar I = 2 mode, and the A p\ for the quadrupolar I = 4 
mode. The 2 f\ dipolar mode is only visible in the spectrum of 
v e . The v° spectrum also shows a peak at 2.45 kHz, which does 
not correspond to any of the known low frequency modes. It is 
possible that its origin might be due to some form of non-linear 
coupling, or perhaps an effect of the free atmosphere. However, 



as already noted in the ID radial case, the presence of a freely 
evolving atmosphere does not seem to affect the frequency of the 
modes, at least within the accuracy of our temporal series. It is 
interesting to recall that inertial modes can be also excited with 
a continuum spectrum extending in a frequency range between 
and 2Q (Font et al. 2002). The peak at 500 Hz, corresponds in 
fact to the rotation frequency of the star. 

We have repeated the same analysis for the model BU8, with 
the same simulation setup. This represents a rapidly rotating 
case, close to the mass shedding limit, and it is a more demand- 
ing test than the mildly rotating case, BU2. Results are shown 
in Fig. 5. Many of the same considerations, as for the previous 
BU2 case, still apply and will not be repeated. Typical devia- 
tions in the value of the density are larger (about a factor 2) than 
in the previous case. There is also evidence for a secular drift (in- 
crease) of the central density, whose average value at the end of 
the run is a factor 10 3 higher than at the beginning. Again this 
drift seems to be twice what is observed in the BU2 case, sug- 
gesting a drop in accuracy at higher rotation rates. We also show 
a Fourier transform of various fluid quantities in time. These are 
all monitored at the same selected location: r — 3.0 and 8 = n/4. 
It is clear that in this case the first I = mode F is by far the 
most strongly excited. However, power in the first I = over- 
tone Hi, at the quadrupole I = 4 mode i p\, and I = 2 overtone 
2 pi, is also present. Little power seems to be present at the I = 2 
2 fx mode. The time evolution of the v e component is quite noisy: 
there is no evidence of power at the H mode frequency, the peak 
at 1 .7 kHz might be a poorly resolved 2 f\ mode, or a contam- 
inated inertial mode close to 2Q frequency. The power that is 
visible around 700 - 800 Hz is likely due to inertial modes, and, 
analogously to the BU2 case, its frequency corresponds to the 
rotation frequency of the star. 

5.5. Stability of differentially rotating magnetized equilibria 

We present here the evolution of an equilibrium configura- 
tion, which is both differentially rotating and contains a strong 
toroidal magnetic field. The conditions for such equilibrium and 
how to build it have been described in Sec. 4. There are no sim- 
ilar tests presented in the literature, and, as a consequence, no 
reference against which to compare our results. However, this 
setup allows us to investigate a strongly magnetized case, where 
we expect the magnetic field to have significant dynamical ef- 
fects. Tests in the literature have, at most, focused on the case of 
weak poloidal fields (Cerda-Duran et al. 2008), for the reasons 
discussed in Sect. 4. However, we are interested here is eval- 
uating a case with a magnetic field that is not simply a small 
perturbation, but that has enough energy to modify the under- 
lying equilibrium. Given that the algorithm described in Sec. 4 
consider only toroidal fields, we limit our analysis to such a case 
and we hope that the proposed test can become a standard, once 
our XNS solver will be made of open use. 

The equilibrium model has the same central density of the 
non-magnetized cases p c = 1.280 x 10~ 3 , and the usual poly- 
tropic law K - 100, n = 1. We adopt a differential rotation pro- 
file with Q. c = 0.02575, A 2 = 70, and a magnetic field charac- 
terized by a magnetic polytropic index m = 1 and coefficient 
K m = 3, corresponding to a maximum magnetic field inside the 
star of =i 5 x 10 17 G, and to a ratio of magnetic energy to total in- 
ternal energy ^0.1 (we name this new model BM). With respect 
to an equilibrium with a very similar rotational structure but no 
magnetic field (model B2: Stergioulas et al. 2004; Dimmelmeier 
et al. 2006) the star is a factor ~ 1.5 larger in equatorial radius. 
The domain, r = [0, 16] and 8 = [Q,n], is covered with 200 
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Fig. 6. Evolution of a stable magnetized solution. The upper panel 
shows a comparison between the initial profiles (solid lines) of the equa- 
torial (green) and axial (blue) densities, the equatorial rotational veloc- 
ity v (magenta) and toroidal magnetic field B (red), and the value of the 
same quantities at at t=1500 (dashed lines). Densities are normalized to 
the initial central value, velocity to its maximum (0.09810) initial value, 
and magnetic field to its maximum initial value too. The insert shows 
the relative difference between the equatorial densities as a function of 
radius. The spike at r ~ 12 is due to diffusive relaxation at the boundary 
between the high density star and low density atmosphere. The lower 
panel shown the variation in time of the central density, at 6 = n/2. 
Note that there is a secular trend of increasing density, with an average 
increase during the simulation run of ~ 1.5 x 10" 3 . 



zones in the radial direction and 100 zones in angle. The star 
is resolved on average over about 120 radial zones. The evolu- 
tion is followed until a final time f max = 1500 corresponding in 
physical units to ^ 7.5 ms. 

Fig. 6 compares the equatorial and axial profiles of the den- 
sity, together with the rotational velocity v = (y^) 1 ^ 2 , and the 
magnetic field B = (B B <Z ') 1/2 , at time t = and t = f max . There 
is evidence for a secular drift (increase) of the central density, 
whose average value at the end of the run is a factor 1.5 10~ 3 
higher than at the beginning, together with typical oscillation of 
similar amplitude. Larger deviations of order of \0 2 character- 
ize the outer layer of the star for r > 10, where the ratio of 
magnetic pressure to gas pressure is larger. This is also visible 
in the broader shear layer that form at the boundary with the at- 
mosphere, compared to non magnetized cases (Figs. 4-5). There 
is also some evidence for a drop in the rotation rate of the core. 
The reason is not clear, it could either be due to some angular 
momentum redistribution from the core to the outer layer of the 
star, or to some convection taking place in the core, possibly ex- 
cited by initial relaxation. It is likely that, due to round-off errors 
and to the initial relaxation, some free energy might be injected 
into the star to power convective motions in regions of marginal 
convective stability as the core. We have however verified that 
the typical poloidal velocities in the core are ~ 10~ 4 , which im- 
plies that it will be necessary to follow the system for a much 
longer time in order to see whether this convective motions are 



stabilized. This also agrees with the results presented by Kiuchi 
et al. (2008), which suggest, even for the case of strong instabil- 
ities, typical convective timescales of a few thousands of M. 

A detailed analysis of the stability of strongly magnetized 
configurations is beyond the scope of this paper. This test how- 
ever demonstrates that stable configurations are preserved for 
many sound-crossing times, in cases involving strong magnetic 
fields too. We have not carried out a full mode investigation, 
as in the previous cases, because there are no known values to 
compare our results to. There is also no reference solution, as in 
the case of RNS for non magnetized rigid rotators. However, al- 
ready from the plot of the central density, it is possible to see an 
oscillatory behavior, with a typical frequency 900 Hz. The am- 
plitude of the oscillations, 10~ 3 , is about one order of magnitude 
smaller that the result of Kiuchi et al. (2008). This frequency can 
be compared with the fundamental mode of a non magnetized 
configuration, with a similar velocity profile (Q. c - 0.02575, and 
A 2 = 70), which is located at 1.37 kHz (Dimmelmeier et al. 
2006). There are also some hints that the first I = overtone 
might have a similarly smaller frequency, even if there is little 
power in it and its identification in the spectrum is not certain. 
A lower frequency for / = modes can easily be understood, 
if one recall that it is well known (e.g. Bucciantini et al. 2003) 
that toroidal magnetic fields behave as gas with adiabatic index 
4/3, for / = perturbations. Compared with the gas, which has 
adiabatic index 2, the presence of strong toroidal magnetic field 
leads to a softening of the generalized EoS of the perturbations, 
corresponding to a lower sound speed, and longer vibrational pe- 
riods. The larger radial extent of the star also contributes, given 
that the period or radial compressive modes scales as the sound 
crossing time. 

5.6. Radial and axisymmetric collapse to a black hole 

The migration test performed in Sect. 5.2 already shows the 
ability of the XCFC algorithm to handle the dynamical evolu- 
tion of very compact configurations. Although that test is failed 
by the original CFC (Dimmelmeier et al. 2002), Marek et al. 
(2006) have shown that it can still be succesfully simulated by 
using a modified version. The superiority of XCFC formulation 
fully manifests itself in cases where the evolution leads to the 
formation of black holes and appparent horizons (AHs) (York 
1989; Baumgarte & Shapiro 2003; Thornburg 2007). To prop- 
erly evaluate the strength of our XCFC solver, we present in 
this section the results of two simulations of the collapse of un- 
stable neutron stars to black holes: a ID collapse of a spheri- 
cally symmetric NS, and a 2D collapse of a rapidly rotating one 
(Bernuzzi & Hilditch 2010; Cordero-Carrion et al. 2009; Baiotti 
et al. 2005). In the ID case the AH is found using a zero-finding 
algorithm, while for the 2D collapse we use a simple minimiza- 
tion algorithm, with the AH parametrized as a surface in terms 
of Spherical Harmonics [for details on the algorithms see Sect 8 
of Thornburg (2007), Sect 6.7 of Alcubierre (2008), and Shibata 
(1997)]. 

For the ID case we consider the same unstable configura- 
tion of the migration test, Sect. 5.2, and following Bernuzzi & 
Hilditch (2010) we add a perturbation of the form 

6 i -0-005 sm(r/r NS ) if r < r NS 

\ if r > r NS v ' 

where tns is the initial radius of the neutron star. The computa- 
tional domain r - [0, 10] is covered by 300 equally spaced radial 
zones. A hot low density atmosphere is set outside the NS, and 
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Fig. 7. Upper panel: the monotonically rising red curves represent the 
evolution of the central density with respect to the intial value pdp c (t = 
0) (as explained in the text the curves are truncated at the formation of 
the apparent horizon), while the monotonically decreasing blue curves 
represent the value of the lapse a at the origin. Solid lines indicate the 
spherically symmetric collapse and dashed lines the collapse of the ro- 
tating NS model D4. Bottom panel: the rising red curves represent the 
radius of the apparent horizon, while the decreasing blue curves are the 
ratio of the rest mass outside the apparent horizon with respect to the 
total rest mass. Again, solid lines indicate the spherically symmetric 
collapse and dashed lines the collapse of the rotating NS model D4. 
The vertical dotted line indicates the time t AH when the apperent hori- 
zon forms. 



let evolve freely. An ideal gas EoS is used, but given that the 
evolution of the collapse does not lead to shocks and dissipative 
heating, results are equivalent to the case of a polytropic EoS. 
Despite our resolution being 300 times worse than the central 
resolution used by Cordero-Carridn et al. (2009), we are able 
to follow the evolution of the system past the formation of an 
apparent horizon. 

In Fig. 7 we show the evolution of a few quantities. The 
apparent horizon forms at a time tAH — 47 to be compared to 
t AH s 48 in Bernuzzi & Hilditch (2010). When the AH forms, 
its location is r = 0.61, the rest mass outside it is 0.28 of 
the initial rest mass, and the value of the lapse at the center is 
a c = 0.028. Our results agree with what is shown in Fig. 3 of 
Cordero-Carridn et al. (2009), and Fig. 8 of Bernuzzi & Hilditch 
(2010). In Fig. 7 we also show the value of the central density. 
However a cautionary remark is here in order about this quantity. 
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Fig. 8. Evolution of the collapse of the rotating unstable equilibrium 
model D4 (density in LoglO units): upper panel, density at t = t AH - 6; 
middle panel density at t = t A H\ lower panel, density at t = t AH + 6. 
The dashed countour indicates the position of the apparent horizon. At 
the end of the simulation a left over disk has formed in the equatorial 
region while the neutron star has been completely accreted in the polar 
region. 



As the collapse proceeds the metric becomes progressively more 
curved in the origin and the density correspondingly steeper, to 
the point that interpolation and round off errors become more 
and more relevant. In Cordero-Carrion et al. (2009), despite their 
much higher resolution, already at t ~ tAH ~ 10 noise is present in 
the plot of the central density, which is no longer monotonically 
increasing. In our case we find that, as the system approaches the 
formation of the AH, the confidence and accuracy with which we 
can derive the central density at r — 0, by extrapolating the val- 
ues on the numerical grid, rapidly decreases. At the time the AH 
forms we have evaluated that the extrapolation accuracy is such 
that the possible error on the density is ~ 10%, and it rapidly 
increases afterward. For this reason, we truncate the central den- 
sity plot at the formation of the AH. It should be reminded that 
the standard practice to handle systems after the formation of 
an AH is to excise the region inside the AH itself, excluding it 
from the computational domain (Baiotti et al. 2005; Hawke et al. 
2005). The precise value of the cental density is obviously not an 
important quantity, as far as the global evolution of the system is 
concerned. Indeed, quantities that depend on global properties of 
the system like the AH radius, Tahi the rest mass left outside it, 
the value of the lapse at the center (which depends on the global 
distribution of matter in the domain), all agree very well with 
what has been previously found in the literature. 

For the 2D case we consider the model D4 of Baiotti et al. 
(2005) and Cordero-Carrion et al. (2009). This corresponds to 
a uniformly rotating neutron star with a central density p c = 
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Fig. 9. Evolution of the toy collapse model. Panels show the rest mass density and the mangetic field lines (represented by arrows) at three different 
instants of the evolution. Left panel: t = 50, initial aspherical collapse. Note that the star is collapsing along the axis but expanding at the equator. 
Middle panel: t = 100, bounce. The star has turned into a disk (aspect ration ~ 1/10). Right panel: t = 200, later expansion. The density is now 
about 15% of the initial value. A density bump is formed distorting the magnetic field. 



3.116 x 10~ 3 , a rotation rate Q = 0.0395, a gravitational mass 
M = 1.86, an equatorial radius r e = 7.6, and an ellipticity 
r p /r e = 0.65. Following Cordero-Carrion et al. (2009) the col- 
lapse is triggered by reducing the pressure 2% with respect to 
the equilibrium value. The computational domain r = [0, 10], 
6 = [0, 7r] is covered by 200 equally spaced radial zones, and 100 
equally spaced angular zones. As usual, a hot low density atmo- 
sphere is set outside the NS, and let evolve freely and an ideal 
gas EoS is used. Again, despite our radial resolution being 50 
times worse than the central resolution used by Cordero-Carrion 
et al. (2009), we are able to follow the evolution of the system 
past the formation of an AH. 

In Fig. 7 we show the evolution of a few quantities. The 
AH forms at a time t^H - 126, to be compared to tAH —130 
in Cordero-Carrion et al. (2009), its location at that moment is 
r = 0.75, the rest mass outside is 23% of the initial rest mass, 
and the value of the lapse at the center is a c = 0.025. Our re- 
sults agree with what is shown in Fig. 3 of Cordero-Carrion 
et al. (2009). The same considerations stated above for the cen- 
tral density still apply. Interestingly, we are able to follow the 
evolution of the model D4 for a much longer time after the for- 
mation of the AH with respect to Cordero-Carrion et al. (2009). 
Fig. 8 shows the evolution of the NS as the AH forms and grows. 
The middle panel represents the system at the moment of the for- 
mation of the AH, as in their Fig. 4 (warning: km units were used 
on the axes). The lower panel shows the density at t — t^H + 6. 
Clearly a disk has formed as in Baiotti et al. (2005): the NS has 
been completely accreted inside the AH in the polar region up 
to a latitude of around 30° while matter is still present in the 
equatorial region up to r ~ 3. 

5.7. Toy collapse 

Core-collapse simulations are often presented as a test for code 
performances. However, they often aim at simulating realistic 
systems and involve complex physics: to the complexity of ex- 
act MHD solutions in a dynamical metric (Bocquet et al. 1995), 
they usually add sophisticated initial conditions from stellar evo- 
lution, tabulated EoS, and neutrino transport. While all these el- 
ements are undoubtedly important in the study of the physics of 
core collapse, it might be questioned wether they are truly neces- 
sary to evaluate the performances of the metric solver. Moreover, 



the diversity in setup, EoS, and transport algorithms makes a di- 
rect comparison among different codes quite difficult. Such sim- 
ulations are also hard to reproduce, if, for example, initial con- 
ditions or tabulated EoS are not easily and freely available, or 
require specific code implementations. Moreover, results them- 
selves often show the growth of turbulence (convection) and 
MRI (Cerda-Duran et al. 2008), for which comparison should 
be done in phase space, and not on selected snapshots at arbi- 
trary times. In an attempt to construct a simple run with fully de- 
fined initial conditions, a simple EoS, and no complex physics, 
we have designed a novel test run of what we call a toy collapse. 
However, this incorporates some of the important elements of a 
fully 2D non-linear evolution typical of a more realistic collapse 
scenario, and it also includes a poloidal magnetic field. 
The chosen initial conditions are the following 

p = io- 4 

p = max[7t/3 (10 2 - r 2 )p 2 , 10~ 9 ] if r < 10 (95) 
= a' 1 (0.025+/?*) 

p = 10~ 7 

p = 10- 9 ifr>10 (96) 
= 0.0 

corresponding to a rigidly rotating uniform sphere. Note that the 
kinetic pressure is 1/2 of the hydrostatic value, for the same 
matter distribution, in Newtonian gravity and with no rotation. 
To this configuration a poloidal magnetic field is added using a 
generator potential 

A,p = max[3 - >/[8(r - 6) 2 + 36(0 - tt/2) 2 ], 0.], (97) 
yielding the components 

& := y" 2 B r = d g A^ (98) 

S e := y 1/2 B B = d r A^, (99) 

where we recall that y 1 ^ 2 = iff 6 r 2 sin 6. 

The initial conditions are found by solving the XCFC equa- 
tions together with the above conditions for the fluid and mag- 
netic variables. This provides us with a self-consistent set of con- 
served variables and metric coefficients, corresponding to our 
choice of primitive variables (note that our definition of velocity 
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implicitly involves the metric). The evolution of this configu- 
ration incorporates various elements of strongly dynamical and 
aspherical collapse. The original uniform sphere collapses pref- 
erentially along the polar axis, where there is no centrifugal sup- 
port, while it tends to expand at the equator, where rotation is 
stronger, forming a strongly oblate disk. Due to the stiff EoS, 
the system bounces when the central density has increased of 
about 50%, then the structure re-inflates, driving a shock into 
the lower density of the surrounding atmosphere. During the col- 
lapse the poloidal field is stretched, compressed, and twisted by 
rotation, but remains confined inside the star. In fig. 9 we show 
three snapshots of the evolution, that we describe as follows: the 
quasi-spherical initial collapse, the bounce, and the final struc- 
ture. The numerical domain, r — [0, 35] and 9 = [0, n], is made 
up by 500 zones in the radial direction and by 150 zones in an- 
gle. The initial high density sphere is resolved over about 150 
radial zones. The evolution is followed for a time t mdx = 200 
corresponding in physical units to ^ 1 ms. 

6. Conclusions 

In this paper we have upgraded the Eulerian conservative high- 
order code for GRMHD (ECHO: Del Zanna et al. 2007) to 
dynamical spacetimes. We have chosen a fully constrained 
method and conformal flatness for the Einstein equations, and 
in particular we have built a numerical solver based on the 
extended conformally flat condition scheme (XCFC: Cordero- 
Carrion et al. 2009) for the elliptic equations providing the 
metric terms. This is known to improve on the previous CFC 
formulation (e.g. Wilson & Mathews 2003), both because of 
the hierarchical nature of the equations to solve (the ellip- 
tic equations are fully decoupled), and because local unique- 
ness of the solution is ensured even for highly dynamical 
non-linear cases. Our novel scheme is here named X-ECHO 
and we also present a code to produce self-consistent initial 
data (metric terms and GRMHD quantities) for polytropic, dif- 
ferentially rotating, relativistic (neutron) stars with a toroidal 
magnetic field, named XNS. This code is publicly available 
at http : //sites . google . com/site/niccolobucciantini/xns 
and we hope it will provide useful benchmark initial data for 
NS evolution or core collapse in the magnetized case. 

Both X-ECHO and XNS work in spherical coordinates of 
the conformally flat metric and axisymmetry is assumed. The 
2D metric solver for the elliptic equations (Poisson-like scalar 
PDEs and vector Poisson PDEs) use a mixed technique: spectral 
decomposition in spherical harmonics (or vector spherical har- 
monics) in the angular direction and finite-differences leading to 
the inversion of band-diagonal matrices in the radial direction. 
This is achieved on the same numerical grid used for evolv- 
ing the fluid/MHD quantities, thus avoiding the need for inter- 
polation over different meshes. We fully test the codes against 
known problems involving fluid configurations in dynamical 
spacetimes, basically ID and 2D evolution and vibration modes 
for NS configurations, migration to stable branches, ID and 2D 
collapse of an unstable NS towards a BH, including the forma- 
tion of an apparent horizon, and we propose a couple of novel 
GRMHD test problems, the evolution of a differentially rotat- 
ing magnetized NS and a toy collapse simulation in the pres- 
ence of poloidal magnetic fields. The metric solver is fast, and, 
on the cases we have tested, the CPU time required to solve 
the XCFC system, is comparable with the time taken to up- 
date the MHD fluid quantities, despite the use of fast Riemann 
solver (HLL/HLLC). For higher than second order reconstruc- 
tion tecniques, the computational time is always dominated by 



the HD/MHD module. The code has been validated against pre- 
vious results obtained with both free-evolution and fully cos- 
trained schemes, and with the linear theory for perturbations. 
Performances in the presence of strong magnetic fields, violently 
dynamical configurations, large deviation from sphericity, and 
even apparent horizons, show that the method and its implemen- 
tation with the HD/MHD module are stable in the situations of 
interest. Moreover we have shown that, in the case of NSs sur- 
rounded by a low density atmosphere, there is no need to apply 
any reset procedure to the atmosphere itself, which can be let 
evolve freely, without altering the stability or the results of the 
simulations. 

For the immediate future we plan to investigate the stability 
and to find the characteristic vibrational modes of a set of mag- 
netized NS configurations, with and without differential rotation, 
both for stable and unstable megnetic profiles. We also plan to 
investigate the growth of magnetic field due to MRI in differen- 
tially rotating NS, which might be of some interest to explain 
the late flaring activity that is observed in long duration GRBs, 
and that, within the millisecond-magnetar model, is commonly 
attributed to bursty magnetic activity in the cooled and convec- 
tively stable NS. However, the final goal is to include a more 
realistic treatment of the microphysics, going beyond the simple 
ideal gas law implemented here for reproducibility of the nu- 
merical tests, especially as far as neutrino heating is concerned, 
and possibly to couple our code with a transport algorithm for 
neutrinos as required for collapse calculations. This will allow 
us to study in details the magnetized core-collapse scenario, and 
to investigate the role of a strong magnetic field in shaping and 
regulating the collapse. More in particular, we also plan to derive 
from magnetized collapse simulations a more realistic setup for 
the (long) GRB model recently proposed by Bucciantini et al. 
(2009), where a newly born millisecond proto-magnetar drives 
a GRMHD wind, that, due to confinement of the external stel- 
lar envelopes, collimates relativistic jets escaping the progenitor 
along the poles. 
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